Cahn-Hilliard theory for unstable granular flows 



T.P.C. van Noije and M.H. Ernst 
Instituut voor Theoretische Fysica, Universiteit Utrecht, Postbus 80006, 3508 TA Utrecht, The Netherlands 

A Cahn-Hilliard-type theory for hydrodynamic fluctuations is proposed that gives a quantitative 
description of the slowly evolving spatial correlations and structures in density and flow fields in the 
early stages of evolution of freely cooling granular fluids. Two mechanisms for pattern selection and 
structure formation are identified: unstable modes leading to density clustering (compare spinodal 
decomposition), and selective noise reduction (compare peneplanation in structural geology) leading 
to vortex structures. As time increases, the structure factor for the density field develops a maximum, 
which shifts to smaller wave numbers. This corresponds to an approximately diffusively growing 
length scale for density clusters. The spatial velocity correlations exhibit algebraic decay ~ r~ d on 
intermediate length scales. The theoretical predictions for spatial correlation functions and structure 
factors agree well with molecular dynamics simulations of a system of inelastic hard disks. 



PACS numbers: 45.70.-n, 45.70.Qj, 47.20.-k, 05. 
keywords: granular fluids, pattern formation, 
equations 

I. INTRODUCTION 

Granular fluids in rapid flow stand out as an interesting 
and complex many-body problem, on which an impres- 
sive amount of experimental, simulation and theoretical 
results have been gathered in recent years. Advanced 
high speed measuring techniques, used in real experi- 
ments, as well as sophisticated visualization methods, 
used in computer experiments, have produced a lot of 
high quality data, which have led to intuitive and qualita- 
tive explanations. However, the experimental conditions 
are either too complex or the relevant physical parame- 
ters unknown to allow quantitative theoretical modeling 
and analysis. 

On the other hand there exists a large body of analytic 
results obtained from hydrodynamic equations, from ki- 
netic theory and from simplified mathematical models 
(on instabilities in flows, on clustering, on transport 
properties and non-Gaussian features in velocity distri- 
butions) that are not (yet) accessible to measurements in 
laboratory or computer experiments. There are a few no- 
table exceptions, which will be discussed below in so far 
as they are relevant for the subject of the present paper. 

Recently we have proposed a new mesoscopic theory, 
based on fluctuating hydrodynamics with unstable modes 
(Cahn-Hilliard-type theory) that provides detailed pre- 
dictions on spatial correlation functions and structure 
factors, which permit a quantitative comparison with 
molecular dynamics (MD) simulations on fluids of in- 
elastic hard spheres. Preliminary accounts of the results 
were published in Refs. This paper || gives a com- 

prehensive account of the underlying theoretical ideas, 
quantitative predictions and the detailed calculations. A 
subsequent paper g gives a comprehensive account of 
the quantitative comparisons with MD simulations, qual- 
itative explanations where detailed theory is lacking, as 
well as discussions on ranges of validity and unresolved 
issues. 

The inelasticity of the collisions between grains makes 
driven and undriven granular fluids behave very differ- 
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ently from atomic or molecular fluids. A dramatic differ- 
ence with elastic fluids is that granular fluids lose kinetic 
energy through inelastic collisions and cool if energy is 
not supplied externally. In a thermodynamic sense, gran- 
ular fluids should be considered as 'open' systems with 
an energy sink, created by the inelastic collisions. This 
collisional dissipation mechanism introduces several new 
time and length scales, which are often related to insta- 
bilities. In this paper we will focus on spatial correlation 
functions and structure factors, and on the underlying 
instabilities in freely evolving undriven granular fluids. 
There are two different mechanisms for pattern forma- 
tion, one comparable to spinodal decomposition j^], a 
second one comparable to peneplanation in structural 
geology ||. Moreover, we will discuss the importance 
of different intermediate scales, related to viscosity, heat 
conductivity and compressibility, which are controlled by 
the inelasticity. 

The instabilities in undriven granular fluids have been 
studied by several authors |0-||,0,|Jl3[Q . Goldhirsch 
and Zanetti were the first to perform molecular dy- 
namics (MD) simulations of an undriven two-dimensional 
(2D) system of smooth inelastic hard disks and observed 
the spontaneous formation of density clusters. The sys- 
tem is unstable against spatial density fluctuations, so in- 
homogeneities in the density field (clusters) slowly grow 
to macroscopic size. However, before this happens, the 
granular fluid, prepared in a spatially homogeneous state, 
remains in a spatially homogeneous cooling state (HCS) 
with a slowly decreasing temperature. Gradually spatial 
inhomogeneities appear in the flow field (vortex patterns), 
and only much later density clusters are being observed. 
In Fig. [j] we show typical snapshots of the momentum 
field and the density field, as obtained in MD simula- 
tions of a system of inelastic hard disks 0^]. The flow 
field develops a large vorticity component and evolves 
into a 'dense fluid of closely packed vortex structures', 
which is still homogeneous on scales large compared to 
the average vortex diameter L v , provided the systems 
are sufficiently large. More detailed information on the 
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dynamics of clusters is available from very recent MD 
simulations EE]. 

The linear stability of undriven granular fluids was 
first analyzed by McNamara || on the basis of hydrody- 
namic equations for granular fluids. In Sec. [o| a detailed 
stability analysis is presented, which will be needed in 
analytic calculations of structure factors and correlation 
functions. The patterns that develop in the flow field cor- 
respond in fact to a relative instability: as momentum is 
conserved in collisions, perturbations in the flow field of 
large enough wavelength decay slower than the typical 
'thermal' velocity, which decays as a result of the inelas- 
ticity of the collisions. In other words, as the system 
evolves, the noise from the thermal motion is reduced, 
which makes the initial long wavelength fluctuations in 
the flow field observable. However, their magnitude re- 
mains bounded by their values in the initial homogeneous 
state. In d-dimensional granular fluids, the modes corre- 
sponding to this relative instability are the (d — 1) trans- 
verse velocity or shear modes, as well as the long wave- 
length longitudinal velocity modes. In Sees. || and III 
we will show how the velocity fluctuations, at finite wave 
numbers, couple to density fluctuations, resulting in a lin- 
ear instability of the density field. Unlike the flow field 
instability, the density instability is an absolute instabil- 
ity, giving rise to macroscopic clustering. 

Comparison of Figs. [I] and || shows how the spatial in- 
homogeneities in density and flow field, slow down the 
cooling. Initially the system follows Haff's cooling law 
[p~6|| , as will be briefly discussed in the next section. Once 
the inhomogeneities have become important, the homo- 
geneous cooling law breaks down at a time r c , estimated 
in Fig. ^ and in Ref. |l4|, and the system crosses over 
to the nonlinear clustering regime with a slower energy 
decay. This happens the more rapidly, the larger the 
inelasticity e. In Ref. 0, a mode coupling theory is de- 
scribed, which perturbatively takes the effect of the in- 
homogeneities on the energy decay into account and em- 
ploys the decay of long wavelength fluctuations. This the- 
ory gives quantitative explanations of the energy decay 
observed in Fig. || at short and long times provided that 
the inelasticity (a > 0.6) and packing fraction (</> < 0.4 
in 2D) are not too large. A more general mode coupling 
theory is developed in Ref. p7[ . 

In order to understand and quantify the patterns of 
Fig. |l|, we study the time dependent spatial correlation 
functions of the corresponding fluctuations, using Landau 
and Lifshitz's theory of fluctuating hydrodynamics p8[ , 
adapted to dissipative hard sphere fluids, i.e. adapted to 
the presence of unstable modes. We present a theory 
that describes the buildup of spatial correlations in the 
density and flow field in the initial regime, where the in- 
homogeneities are governed by linear hydrodynamics, i.e. 
linearized around the homogeneous cooling state (HCS). 
To motivate and explain the theory we compare a freely 
evolving granular fluid with spinodal decomposition pi, 
where the observed phenomena are similar in several re- 
spects. 



In spinodal decomposition a mixture is prepared in a 
spatially homogeneous state (analogous to the homoge- 
neous cooling state) by a deep temperature quench into 
the unstable region. In this region long wavelength com- 
position fluctuations are unstable (analogous to density 
fluctuations in granular fluids) and lead to phase separa- 
tion and pattern formation (analogous to the formation 
of density clusters and vortex patterns). This instabil- 
ity and the concomitant pattern formation are explained 
by the Cahn-Hilliard (CH) theory |^|, where the time 
dependent structure factor for composition fluctuations, 
S nn (k, t), is calculated on the basis of a macroscopic dif- 
fusion equation. In the CH theory the unstable composi- 
tion fluctuations are described by a fc-dependent diffusion 
coefficient, which is negative below a certain threshold 
wave number. 

The Cahn-Hilliard theory has also been used in the the- 
ory of two-dimensional turbulence (see Frisch jl9| and 
references therein), where the behavior of the fluctua- 
tions in the flow field of incompressible fluids has been 
described in terms of negative eddy viscosities. Vorticity 
modes with negative effective viscosities have also been 
used by Rothman to study vortex formation in lattice 
gas cellular automata pOfl . 

The instability of undriven granular fluids also differs 
in many details from spinodal decomposition. The for- 
mer is a slow, the latter a fast process. Consequently, 
the Cahn-Hilliard theory in spinodal decomposition only 
describes the onset length and time scales of phase sep- 
aration. As the formation of vortices and clusters in 
undriven granular fluids is a rather slow process, the 
present theory is expected to give a good description up 
to times which are rather large (see Figs, [l] and |^) on 
time scales that can be reached in molecular dynamics 
simulations, provided the inelasticity and the density are 
not too large. 

The dynamics of fluctuations in granular fluids, say, in 
density, 5n(r,t) — n(r,t) — (n), a nd flow field, u(r,i), 
have hardly been studied ||7 yl2|Jl3|| , in sharp contrast to 
the average behavior, as witnessed by a large number 
of publications. We first recall that fluctuations are ab- 
sent in standard hydrodynamic as well as in Boltzmann- 
Enskog-typc kinetic equations which are based on molec- 
ular chaos (mean field assumption). 

The objects of interest in this article are the equal- 
time spatial correlation functions of the fluctuating hy- 
drodynamic fields 5a(r,t) = a(r,t) — (a), where a,b — 
{n,T,u a } with a = x,y,... denoting Cartesian compo- 
nents, 



G ab (r,t) = V~ 1 / dr'{5a(r + r',t)Sb(r',t)). 



(1) 

The structure factors are the corresponding Fourier 
transforms, 

S ab (k, t) = J dr exp(-zk • r)G a6 (r, t) 

= V- 1 (Sa(k,t)Sb{- 1 k,t)), (2) 
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where Sa(k,t) is the Fourier transform of 6a(r,t), and 
V = L d is the volume of the system. Goldhirsch et 
al. JtJ initiated MD studies of 5 ran (k, t) and S pp (k,t) = 
Y] S aa (k, t), and related in a qualitative way the struc- 
ture at small k to the most unstable shear modes. They 
presented a nonlinear analysis to explain the enslaving 
of density fluctuations by the vorticity field. This analy- 
sis reveals that the length scale associated with the late 
stages of nonlinear clustering is of the order £j_ ~ l / y/e, 
where e = 1 — a 2 measures the inelasticity and Iq is the 
mean free path. Brey et al. have also studied the nonlin- 
ear response of the density field to an initially excited k 
mode in the transverse flow field [^lj . 

The most important function to describe the cluster- 
ing instability is the structure factor S nn (k,t). A first 
step in its theoretical understanding has been given by 
Deltour and Barrat fj^] . These authors have shown how 
its growth rate in the linear regime is determined by the 
most unstable long wavelength part of the heat mode, in 
which the density couples to longitudinal velocity pertur- 
bations. 

The main goal of the present paper is to study the ve- 
locity and density correlation functions in undriven gran- 
ular flows and to demonstrate the importance of internal 
fluctuations. In Sees. [II and IV wc discuss two different 



mechanisms for pattern formation, one leading to den- 
sity clusters (spinodal decomposition), and one leading 
to vortex structures (peneplanation) , and we we calcu- 
late the structure factors and show that a mesoscopic de- 
scription in terms of fluctuating hydrodynamics Jl8| gives 
quantitative predictions for the spatial correlations func- 
tions over large intermediate time intervals, controlled by 
linearized hydrodynamics. 

In Sec. [v| an analytic description of the correlation 
function G a p(r,t) of the components u a (r,t) of the flow 
field will be given, based on fluctuating hydrodynamics 
and the assumption of incompressible flow. This theory 
was presented in Ref. Q and yields predictions, includ- 
ing long range tails ~ r~ d in d-dimcnsional fluids, that 
for nearly elastic particles (a > 0.9) agree well with two- 
dimensional MD simulations up to large distances. As 
the transverse velocity fluctuations of an incompressible 
fluid do not couple to the density fluctuations in a linear 
theory, this theory gives no information on S nn (k,t). A 
quantitative theory for the magnitude of S nn (k,t) and 
its k dependence in the full range of hydrodynamic wave 
numbers was given in Ref. Q and is described in the 
second part of Sec. [v|, where the assumption of incom- 
pressibility is dropped and the description is extended 
to the general compressible case, based on the full set 
of fluctuating hydrodynamic equations. The results are 
also compared with com put er simulations. We end with 
some conclusions in Sec. [V|. 



II. HYDRODYNAMIC STABILITY 

In the present section we discuss the macroscopic hy- 
drodynamic equations for inelastic hard sphere (IHS) flu- 
ids. Fluctuations will be considered in the next section. 
We assume that IHS hydrodynamics for weakly inelastic 
systems can be described by the standard hydrodynamic 
conservation equations supplemented by a sink term T in 
the temperature balance equation, 

d t n + V • (rai) = 



9 ( u + U'Vu= -V • n 

P 

d t T + u-VT = -- ^(V J + II : Vu) - T. 

an 
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Here p = mn, u the flow velocity, and ^dnT the ki- 
netic energy density in the local rest frame of the IHS 
fluid. The pressure tensor n„^ = p5 a p + SH a/ 3 contains 
the local pressure p and the dissipative momentum flux 
STl a j3, which is proportional to V Q u^ and contains the 
kinematic and longitudinal viscosities v and vi, defined 
below Eq. (Af ) of App. |A|. The constitutive relation for 
the heat flux, J = -kVT, defines the heat conductiv- 
ity k. The above equations can be justified to lowest 
order in the inelasticity e = 1 — a 2 p^ , p3[ , and will be 
used for weakly inelastic systems, where the pressure is 
assumed to be given by its value for elast ic h ard spheres, 
p = nT[l + 2b{n)] with b(n) given below ( |Al| ) in App. [a|. 
Note that 0(e) contributions are being neglected here, 
because the Enskog-Boltzmann equation for IHS fluids 
predicts p = nT[l + (1 + a)b(n)]. Similarly, the trans- 
port coefficients v, v%, and k are assumed to be given by 
the Enskog theory for a dense gas of elastic hard spheres 
(EHS) or disks p4| . In the present paper, corrections of 
0(e) to the leading nonvanishing behavior of thermody- 
namic properties and transport fluxes will be neglected. 
More general expressions for transport coefficients have 
been derived in Refs. p5|-p9|]. They contain an addi- 
tional contribution, — /iVn, in the heat flux, where the 
new transport coefficient jj, ~ e vanishes for small inelas- 
ticity. 

Kinetic theory provides an exact expression for the col- 
lisional dissipation rate T, which represents the average 
energy loss through inelastic collisions. It can be derived 
from the microscopic energy loss per collision. An ex- 
plicit derivation can be found in Refs. [ fl6|]25|]28[ ] . For the 
present purpose however, a phenomenological derivation 
suffices, which proceeds as follows. On average, a particle 
loses per collision an amount ~ eT of its kinetic energy 
and per unit time an amount ~ ewT, where u> is the aver- 
age collision frequency p4| . This argument gives (apart 
from a numerical factor l/d) 



2-fouT, 



(4) 



where 70 = (1 — a 2 )/2d. In the present theory the colli- 
sion frequency is assumed to be given by Enskog's theory 
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for dense hard sphere fluids |24|] , and quoted explicitly in 
Eq. (|A|) of App. |A|. It is proportional to \/T. 

For an understanding of what follows two properties 
of undriven granular fluids are important: (i) the exis- 
tence of a homogeneous cooling state (HCS) and (ii) its 
linear instability against spatial fluctuations. We first ob- 
serve that the hydrodynamic equations for an IHS fluid, 
initialized in a homogeneous state with temperature To, 
admit an HCS solution with a homogeneous density n, a 
flow field which can be set to zero everywhere, and a ho- 
mogeneous temperature T(i), determined by dtT = —V. 
To solve this equation it is convenient to change to a 
new time variable, defined as dr = u>(T(t))dt, yielding 
T(t) = Tq exp(— 2-for). To find the relation between the 
'internal' time r, which measures the average number of 
collisions suffered per particle within a time t, and the 
'external' time t, we integrate the relation for dr using 
lu ~ VT, with the result 



t = — ln[l + jot /t ] 
7o 



(5) 



In the elastic limit (e — > 0), it is proportional to the ex- 
ternal time, r = t/to, measured in units of the mean 
free time to — 1/lu(T ) at the initial temperature. The 
theoretical prediction (||) agrees well with the computer 
simulations Q until the system crosses over at r c into the 
nonlinear clustering regime, to be discussed below. The 
initial slope of r(t) corresponds to the collision frequency 
a; (To) in the equilibrium state at r = 0. The relation (||) 
introduces a new intermediate (mesoscopic) time scale, 
the homogeneous cooling time t e — to/ Jo- Combination 
of these results yields the slow decay of the temperature, 



T(t) = T exp(-2 70 T) = T /(l 



7o*Ao) 2 , 



(6) 



which is Haff's well-known homogeneous cooling law p6[ . 
There exist extensive verifications of the validity of Haff's 
law in the HCS (r <C r c in Fig. ^) using MD simulations 
d,|0]|l|]. We show below that this HCS solution is lin- 
early unstable once the linear extent L of the system 
exceeds some dynamic correlation length ~ lo/\/t j^Di 
where lo is the (time independent) mean free path. It is 
given by lo — Vo/u>, where vq = \j2T/m is the thermal 
velocity. 

In general the HCS is highly nontrivial, as it exhibits 
correlations between the velocities and positions of differ- 
ent particles. In the 'lowest order' description (for more 
refined approximations see ]|a j2q , B(| ) the HCS corre- 
sponds to an equilibrium state, which is cooling adiabati- 
cally, i.e. with time dependent temperature @). Here ve- 
locity correlations between different particles are absent, 
and position correlations are taken into acco unt through 
the pair correlation function at contact (see (A2) of App. 

0). 

In the present paper, we are interested in the buildup 
of correlations between spatial fluctuations in a sys- 
tem that is prepared in a homogeneous state at an ini- 
tial temperature Tq. It reaches the HCS within a few 



mean free times to- Therefore, we can linearize Eqs. 
(0) around the homogeneous density n and temperature 
T(t) = T /[l + jot /t } 2 , and the vanishing flow field of 
the slowly evolving HCS. The resul ting set of linearized 
hydrodynamic equations, given in (Al) of App. |A|, con- 
tains the Enskog hard sphere transport coefficients, v, 
vi and k, which are proportional to y/T(t), and depend 
therefore explicitly on time. It is again convenient to 
make the transformation dr = u>dt, and introduce the 
dimensionless variables Sn(r, r) = 5n(r,t)/n, u(r, r) = 
u(r,t)/v (t), and ST(r,r) = 6T(r,t)/T(t). In these new 
variables the equations of change for the macroscopic 
Fourier modes Sh(k, r), u(k, t), and ST(k, t), defined 
through 5a(k, t) = Jdrexp(— ik • r)Sa(r, r), become or- 
dinary differential equations with time independent co- 
efficients (valid for klo i5 1), which can be analyzed in 
terms of eigenvalues and eigenvectors, 



dSh 

~dr 

duj_ 

du. 
7fr 



—ikloui 

7o(i-fc 2 ei)«± 



2nT 



ST 



iklo 



1 



2nT\T 



dST 



~l^ + k 2 e T W~iklo(^^j ui 
-j g(n)Sh. 



(7) 



Here we have introduced the time independent correla- 
tion lengths £xj 6 an d Ct, defined by £± = v/coja, 
£f = vijiojo and = 2k/ dnojjo, the isothermal com- 
pressibility xt = (dn/dp) T /n, and for the function 
g(n) = 2[1 + (n/x)dx/dn] we refer to Eqs. ( |Al| ) and 
( |A2j ) of App. |A| The subscript _L in the equation for 
u± refers to any of the (d — 1) directions perpendicular 
to k, and the subscript I denotes the longitudinal direc- 
tion along k. The validity of the hydrodynamic Eqs. (0) 
is restricted to wave numbers k -C 2-k/Iq to guarantee 
separation of kinetic and hydrodynamic scales, and to 
k <C 2tt/(t, where a is the diameter of a disk or sphere, 
to guarantee that the Euler equations involve only local 
hydrodynamics. Therefore, the validity of the hydrody- 
namic Eqs. (0) are restricted to long wavelengths, satis- 
fying k <C min{27r/Zo, 27r/cr}. In matrix representation 
we write the above equations as 



JUa(k,T)=M(k)<5a(k,r), 



(8) 



where components of a are labeled with {n, T, I, _L}, and 
the hydrodynamic matrix M is defined by Eqs. ([?]). The 
linear stability of (g) was first investigated in Refs. 0,^). 
Its eigenvalues or dispersion relations Ca(£) and corre- 
sponding eigenvectors are analyzed in App. |X| The C,\(k) 
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were first calculated numerically in Ref. M , using trans- 
port coefficients somewhat different from those obtained 
from the Enskog theory. Typical results of our calcu- 
lations are shown in Fig. [| The most striking feature 
is that there are two eigenvalues, £j_ and that are 
positive below the threshold values k* ± and k* Hl i.e. two 
linearly unstable modes with exponential growth rates. 
The spectrum has different characteristics for smaller and 
larger wavelengths. In the dissipative range || (klo <C 70) 
all eigenvalues are real; propagating modes are absent. 
Around kl ~ 0(70), two eigenvalues become complex 
conjugates and the corresponding (sound) modes become 
propagating. In the standard range (klo ^70), compres- 
sion effects and sound propagation, which are O(klo), 
dominate dissipation, whereas heat conduction, which 
is 0(fc 2 Z 2 ), becomes dominant only in the elastic range 



(klo ^ V^o)' where 7o is assumed to be sufficiently small 
so that 70 <C -y/To- In the latter range, the dispersion re- 
lations and eigenmodes resemble those of an elastic fluid. 
In App. |A| the full set of eigenvalues (\(k) and eigen- 
modes is analyzed. Here we only analyze the (d— I)-fold 
degenerate transverse velocity or shear modes, which are 
decoupled from the remaining modes. In fact, Eqs. (Q) 
show already that Cl(&) = 7o(l — k 2 ^), which describes 
an unstable mode for k < k\ — l/£x- 
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Similarly, the heat mode is unstable for k < k H 
1 /£rj as shown in App. |X| In the dissipative range (klo •C 
70), the heat mode is given by Cn(k) — 7o(l — k 2 £, 2 ) where 
£|| is calculated in App. A|. In this range the heat mode 
is a pure longitudinal velocity fluctuation. We point out 
that £j_ diverges as l/\fe for small e, while £|| ~ 1/e [sec 
(A7)]. As a consequence the correlation lengths £j_ and 
£y are well separated for small inelasticity, as shown in 
Fig. I 

Furthermore, we observe that the instability of shear 
and heat mode is a long wavelength instability. As a con- 
sequence, for finite systems effects of the boundaries are 
important, and the various instabilities are suppressed 
in small systems. When using periodic boundary con- 
ditions, the instabilities are suppressed if £; m in = 2ir/L 
is larger than k± or k* H . When decreasing the system 
length L = V 1 ^ at fixed inelasticity, first the heat mode 
will become stable (k* H < k m i n < k±). In this range the 
density (coupled to the heat mode) is linearly stable, and 
density inhomogeneities can only be created via a nonlin- 
ear coupling to the unstable shear mode 0. Decreasing 
the system size even further (k*\_ < fc m i n ) will stabilize the 
shear mode and thus the HCS itself. In the next section 
we present a mesoscopic theory to describe the dynamics 
of the long wavelength fluctuations in the system. 



III. MESOSCOPIC HYDRODYNAMICS 

In this section we study spatial fluctuations, around a 
reference state, of the hydrodynamic fields, which leads to 
a Cahn-Hilliard-type theory M for the structure factors. 



The dynamics of these fluctuations can be described by 
the fluctuating hydrodynamic equations jl8), obtained 
from the nonlinear equations ([|) by adding fluctuation 
terms to the momentum and heat current, denoted by II 
and J respectively. These currents II and J are consid- 
ered as Gaussian white noise, local in space, and their 
correlations are determined by some appropriately for- 
mulated fluctuation-dissipation theorem for the reference 
state. 

In elastic fluids the reference state would be the ther- 
mal equilibrium state. In driven systems it would be 
a nonequilibrium steady state (NESS). In the present 
case the reference state is the slowly evolving homoge- 
neous cooling state. In lowest approximation p8[ | it may 
be considered as an adiabatically changing equilibrium 
state with a constant density, a vanishing flow field and 
a time dependent temperature, described by Haff's law 
@. The basic extension required for application to IHS 
fluids, is the assumption that the fluctuation-dissipation 
theorem also applies to the HCS with an adiabatically 
changing temperature T(t). This assumption relates the 
noise strengths to the transport coefficients through |lj] 

(tl a0 (r, i)n 7 a(r', t')) = 2T[r](5 ai 6ps + <U<W) 

+(C - ^)<WM<5( r - r'W - 
(j a (r,t)j p (r',t')) =2KT 2 5 a/3 6(r-r')8(t-t'). (9) 

For dimensional reasons, the transport coefficients in sys- 
tems with hard sphere type interactions, like IHS, are 
proportional to ^T(t). 

In the present theory for nearly elastic fluids we lin- 
earize the nonlinear Langevin equations, obtained from 
(||), around the HCS. By applying the transformations 
introduced above Eqs. (|7|) we obtain the dynamic equa- 
tions for the rescaled variables, in the form of a set of 
Langevin equations with constant coefficients M(k), i.e. 



d_ 



Sa(k, t) = M(k)<5a(k, r) + f (k, t), 



(10) 



where f represents the rescaled internal fluctuations in 
the momentum and heat flux. 

The equation of motion for the matrix of rescaled 
structure factors, S a b(k, t) = V~ 1 (5a(k,T)5b(—'k,T)}, 
can now be derived from Eq. (|l0|), and yields 

^-S(k, t) = M(k) • S(k, t) + S(k, r) • M T (-k) + C{k), 

OT 

(11) 

where M T is the transpose of M. It is to be solved 
for given initial values S(k, 0). In terms of the rescaled 
variables the Gaussian white noise (^|) has the standard 
from with constant coefficients, given by the covariance 
matrix, C(k), with 

V- 1 (f a (k,r)f b (-k,T'))=C ab (k)S(r-r'). (12) 
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It is diagonal with nonvanishing elements 
C TT = 47 fc 2 £f./dn 

C^i_ = 7ofc 2 el/n, 



(13) 



as follows from (m and the definitions of the correlation 
lengths below ([tJk The formal solution of ([ll]) for the 
matrix of rescaled structure factors is then 

S(k,r) = exp[M(k)r] • S(k,0) • exp[M T (-k)r] 

dr' exp[M(k)r'] ■ C(fc) • exp[M T (-k)r']. (14) 



At the initial time (r = 0) the system is prepared 
in a thermal equilibrium state of elastic hard spheres 
with density n and temperature To. Consequently, all 
elements of S(k,0) are known. Moreover, the evolu- 
tion equations ( |Tl| ) and (0) are only valid for k <C 
{2tt/Iq, 2ir/a}. So the initial values S(k, 0) are only 
needed for ko~ <C 2n, where they are given by their lim- 
iting values 1 as k — ► 0. The nonvanishing 5 a fc(k, 0) with 
a,b = {n, T, I, _L} are given by their equipartition values 
for elastic hard sphere fluids, 



S nn (k,0) = Txt 
2 

dn 



S TT (k,0) 
S u (k,0) = S_ 



1 

2n' 



(15) 



The above equations can be solved numerically or ana- 
lytically. For a theoretical analysis it is more convenient 
to study the deviations from the initial values in thermal 
equilibrium, defined as 



S + (k,r) = S(k,r) - S(k,0). 



(16) 



The reason is that S(k, r) itself at larger k values also ap- 
proaches S(k, 0) by the imposed fluctuation-dissipation 
theorem (|j) or (p^). Inspection of ( |lT| ) for kl 3> ^/To 
(elastic range) shows that all k independent terms in the 
matrix M(k) in (Q) can be neglected and the hydrody- 
namic matrix reduces to the elastic one, E(k). Conse- 
quently, Eq. ([ll]) reduces to 

E(k) • S(k, 0) + S(k, 0) • E T (-k) + C(fc) = 0, (17) 



as can be verified using ( |13| ) and (llq). 

Subtracting this equation from (O) yields an equation 
of a similar form as (O), 



B(fc), 



|-S+(k,T)-M(k).S + (k,r) 

+ S+(k,r) -M T (-k) 

with a source term 

B(fc) = [M(k) -E(k)] • S(k,0) 

+ S(k,0)-[M T (-k)-E T (-k)]. 

Its nonvanishing matrix elements are 

B n T = B Tn = -j g(n)S nn (k,0) 
B TT = -2'y §TT(k,0) 

Bu=2 lQ Su(k,0) 
B ±± =2 7o 5 ±± (fc,0), 



(18) 



(19) 



and the formal solution of 
S+(k,0) = becomes, 



(20) 

with the initial value 



S+(k,r) 



dr' exp[M(k)r'] ■ B(fc) • exp[M T (-k)r']. 

(21) 



The spectral decomposition (A5) of App. ^ allows us to 
write the rescaled structure factors for a,b = {n,T, 2, _L} 
as 



A/; 

exp[(C A + Cm)t] - 1 
Ca + Cm 

where we have introduced 

B x ,(k) = (v A (k)|B(fc)|v„(-k)) 



(22) 



(23) 



and the eigenvalues CaW depend only on |k| (see App. 
^). The results (14) or (22) yield the structure factors, 
which can be calculated numerically or analytically. Ex- 
pression (14) is most convenient for numerical compu- 
tation, whereas (^2|) is more suitable for our theoretical 
analysis in the long wavelength range. The expressions 
above contain exponentially growing factors describing 
the unstable modes, as in the Cahn-Hilliard (CH) the- 
ory, discussed in the introduction. The present theory 
includes in ([Til ) the rapid microscopic fluctuations, in- 
duced by the Langevin noise C(k), and is equivalent to 
the Cahn-Hilliard-Cook theory M. If the Langevin noise 
is neglected by setting C(fc) =0 in (|l4|), we obtain the 
predictions of the noiseless CH theory, 

S(k, r) = exp[M(k)r] • S(k, 0) • exp[M T (-k)r]. (24) 



1 Note that k should not be set equal to 0, because Sh(0, t) — 
when the total number of particles is fixed. 



G 



or in component form 

S ab (k, t) = J~] wx a (k)w IMb (-k)S\ fl (k) exp[(C\ + C0 T ] 



(25) 



where ^^(k) is defined by (|23|) with B(fc) replaced by 
S(k,0). 

The explicit solutions can be studied using the eigen- 
values and eigenfunctions in different ranges of wave 
numbers. For the very long wavelengths in the dissipa- 
tive range (Hq <C 70 = e/2d) this will be done in the next 
section by /c-expansion at fixed e. In general, the dissipa- 
tive range (Hq <C 70) and the standard range (klg 3> 70) 
are accessible to analysis by rescaling the wavelengths as 
k = e 2 k or k = ^fek respectively, and taking the small-e 
limit subsequently |)). 

A theoretical analysis of the formal results, derived in 
this section, will be presented in Sec. IV. Before con- 



will be presented in Sec. 
eluding this section we present the numerical results of 
( |l4| ) and (|2~i| ) respectively with and without Langevin 
noise, where the relevant equations have been solved with 
Mathematica. The values of 5 a t(k, t) 'with noise' are 
plotted as solid lines in Figs. [| ||, and ^ for different 
components (ab); the values 'without noise' are shown as 
dashed lines in Figs. |5|(a) and ^|. 

The qualitative features of both theories are about the 
same, except that the large-A; limit of the results of ( gj ) 
'without noise' at fixed t does not approach the plateau 
values ([l5|), but vanishes. In the long wavelength and 
large time limit the results of both theories approach each 
other, but at finite k there are substantial differences. It 
should of course be kept in mind that the Cahn-Hilliard 
theory without noise was only designed to explain the 
structure and instabilities on the largest wavelengths. 

It is also of interest to observe that the locations of the 
maxima of S nn in Fig. and those of the minima ('dip') 
of Su in Fig. U approximately coincide. It indicates that 
the density instability is closely connected to the dynam- 
ics controlling Su . This turns out to be the heat mode, as 
we will show in the next section. The negative value of 
SnT{k,r) in Fig. shows that density and temperature 
fluctuations are anticorrelated at all wavelengths. 

An illustration of the comparison with the MD simula- 
tions on systems of inelastic hard disks of Refs. is 
shown in Fig. ||(b) for S nn and in Fig. |^(b) for S± = S±± 
and Su = Su. The agreement between theory and simu- 
lations is in general very good, even down to rather large 
inelasticities (a > 0.6). By comparing the simulation 
results for S± and S\\ in Fig. ^(b) with the numerical re- 
sults of our theory with noise (solid lines) and without 
(dashed lines), we have observed that the agreement in 
the former case extends over the full range of k values, 
whereas in the latter case it is restricted to the small- 
k range. Similar conclusions hold when comparing the 
simulation results for S nn (k,t) in Fig. |B|(b) with the nu- 
merical results in Fig. ||(a) with and without Langevin 



noise included. A more comprehensive comparison with 
MD simulations for different densities and different in- 
elasticities will be given in Ref. where also the range 
of validity of the present theory will be explored. 



IV. DEVELOPMENT OF STRUCTURE AND 
INSTABILITIES 

An elastic fluid in thermal equilibrium does not show 
any structure on hydrodynamic length scales (k < 
min{l/Zo, 1/er}). This means that the hydrodynamic 
structure factors S a b(k) are totally flat, independent of fc, 
as can be seen in (|l5|) . The corresponding hydrodynamic 
correlation functions are short ranged, G a h(r, t) ~ S(r), 
on these length scales. Development of structure on 
length scales above the microscopic scales {lo,&} will 
manifest itself in the appearance of one or more maxima 
or peaks in the structure factors. A linear instability will 
manifest itself in a structure factor that grows exponen- 
tially in time. With these concepts in mind, we analyze 
the structure factors in ( [Ti] ) for the IHS fluid, as we want 
to determine which physical excitations are responsible 
for the features observed in the MD simulations and in 
the numerical solutions. 



A. Structures in the velocity field 

We first consider the simplest case of the the transverse 
structure factor S±(k 1 t) — UQ(i)Sj_(fc, t) with (ab) = 
(_L_L). It describes the transverse velocity or vorticity 
fluctuations u± (k, r) , which are decoupled in (0) from the 
remaining Fourier modes, and satisfy a one-component 
Langevin equation, where the matrix M(k) in (^) re- 
duces to a single number C±(k) = 7o(l — £,±k 2 ). The 
structure factor is readily found from (E2J) and yields 



S+(k,t)=v 2 (t) / dr'S ±± (fc)exp[2a(fc)r' 
Jo 

_ T(t) / exp[27 (l-ejfc 2 )T]-I 
mn \ 



1 - e±k 2 



(26) 



where ^mv 2 (t 



= T(t) = To cxp(— 270T). Combination 
of ( |16| ) and (|26|) yields the complete structure factor, 
S±(k,t) = v%(t)[S±(k,0) + S_f(k,T)], which is plotted 
as the upper solid line in Fig. Wb). It does not grow, but 
slowly decays as S±(k,t) ~ (To/mn) exp(— 2^^\k 2 r) at 
the largest wavelengths. It simply represents vorticity 
diffusion on the 'internal' time scale r, with a diffusivity 
7o£^ = v/lu, where the typical length scales of vortices 
grow |HJ like L v (t) ~ 27r£^\/27o' r ~ 2tt^/vt/uj, which is 
independent of the degree of inelasticity. 

Consider the noiseless CH theory in (pjjl), which 
yields directly S±(k,t) — (T /mn) exp(— 27 ^fc 2 r) [up- 
per dashed line in Fig. ^|(b)] . It agrees with the full the- 
ory ( p6| ) [upper solid line in Fig. ||(b)] only in the small-fc 
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limit, where the denominator (1 — Q_k ) in ( p6[ ) can be 
replaced by 1. At finite k, the differences between both 
theories are relatively large. Moreover, the upper dashed 
line in Fig. ||(b) convincingly shows that the noiseless CH 
theory does not agree with the simulation results. 

Next, we consider the longitudinal structure factor 
S\\(k, t) with (ab) = (11). As subsequent analysis will 
show, it is determined — in an almost quantitative 
manner — by a single mode, the heat mode. The 
fastest growing term in ( |22] ) is (A/i) = (HH), as can 
be seen from Fig. ^. This leads in the structure factor 
S\\(k,t) = VQ(t)S\\(k,r) to a slowly decaying contribu- 
tion. All remaining contributions decay at least as fast 
as vo(t) ~ exp(— 7 r). The slowest decay occurs at very 
long wavelengths, i.e. in the dissipative range (fc^o <C 7o)- 
The rele vant coefficient BnH(k) follows from ([||), (|2C~ 
and (A6) in App. [A], and is given by 



B HH (k) =J2vHa(k)t H b(-k)B ab (k) 

ab 



To 



flow field decrease at larger k values and remain bounded 
for all k and t by their initial equipartition value To / mn, 
which is independent of k. This can be rephrased by stat- 
ing that the flow field exhibits only a 'relative' instability. 

The noise reduction is a direct consequence of the mi- 
croscopic inelastic collision rule, which forces the parti- 
cles to move more and more parallel in successive colli- 
sions. It is this 'physical coarse graining' process that se- 
lectively suppresses the shorter wavelength fluctuations 
in the flow field in an ever-increasing range of wave- 
lengths. Consistent with this picture is also the selective 
suppression of the divergence of the flow field uii(k, t), 
which decays at a much faster rate 7o£jjfc 2 than its rota- 
tional part uj_(k, t) that decays with a rate 7o£j_fc 2 • 

So, noise reduction is the pattern selection mechanism, 
responsible for the growing vortex structures observed in 
Fig. [lj An interesting analog from structural geology of 
the peak formed at k = in reciprocal space, is the for- 
mation of Ayers Rock (Mount Uluru) in the center of 
Australia, as a result of peneplanation by selective ero- 



(27) 



Th e re quired component of the right eigenvector, given 
in (A6), is u)#;(k) ~ 1. Inserting these data in ( f22[) and 
combining it with (|l^) yields for Mq <C 70, 



inn 



has 



where the dispersion relation Cfj(fc) = 7o(l — sii 
been used. This long wavelength approximation for 6*11 
has the same form as the exact result (26) with £j_ 
replaced by In Fig. ||(a) we compare the result 

( p8| ) (dot-dashed lines) with the numerical solution (solid 
lines), presented in Sec. |lll| . The simple analytic small-fc 
approximation (^) captures the global features at small 
k and the plateaux at larger k values quantitatively, but 
is missing the little dip at intermediate k values. 

The behavior of the longitudinal structure factor on 
the largest length scales follows again from Eqs. ( pq ) 
as S\\(k,t) ~ (To/mn) exp(— 270^ fc 2 r). This implies 
that the heat mode on the largest length scales is a 
purely diffusive mode with a diffusivity 7o£f; which is 
much larger than the diffusivity 7o£ 2 _ for the vortic- 
ity (see Fig. ||), and the associated length scale grows 
like Ln(t) ~ 27r^|| v / 2"7q"t. Inspection of the eigenmodcs 
{\v H (k), v#(k)} in ( |A6| ) of App. [A] for k — > shows that 



this diffusive mode is a pure longitudinal velocity field 
£t;(k, t). Its diffusivity 7o£ 2 , given in (JA7|) of App. 
depends for small inelasticities (70 — > 0) mainly on ther- 
modynamic variables, like compressibility and pressure, 
and only slightly on transport coefficients. 

The physical implications of Fig. |^ are quite interest- 
ing. It shows the phenomenon of noise reduction |l4| ] 
at small wavelengths. With increasing time the noise 
strengths Su(k,t) and S±(k,t) of the fluctuations in the 



B. Unstable density structures 

So far, we have seen that the velocity structure fac- 
tors, S±(k, t) and S\\(k, t), develop structure, but remain 
bounded by their initial values. Next we will focus on 
the unstable structure factor S nn (k 1 t), which describes 
the density clustering in undriven IHS fluids. In the com- 
parable case of spinodal decomposition the phase sepa- 
ration is driven by a single unstable mode, the composi- 
tion fluctuations, described by the macroscopic equation 
dtSn(k., t) = zn(k)5n(\i, t), where the growth rate has 
the typical form zu(k) = Afc 2 (l — ^£,%k 2 ). The cor- 
responding structure factor (|^) in the noiseless Cahn- 
Hilliard theory has the form S nn (k,t) ~ exp[2zi)(k)t], 
and exhibits a maximum growth rate at fc max = 1/£d, 
where zu{k) has a maximum. This time independent 
length scale fails to describe the growing length scales of 
the patterns observed in spinodal decomposition 0. 

However, the present theory does explain the growing 
correlation length L c \(t) for the clusters in undriven IHS 
fluids in the early stages of cluster formation, as the nu- 
merical results for S nn (k,t) in Fig. || demonstrate. Its 
maximum at fc max (i) shifts to smaller k values. 

We first consider a naive version of the noiseless Cahn- 
Hilliard theory, proposed by Deltour and Barrat Jl2[ . 
These authors assume that the structure factor S nn can 
be described by the unstable density field, i.e. S nn (k, t) ~ 
S n n(k, 0) exp[2£ff (fc)r], with the growth rate Cn(k) of the 
unstable heat mode. As Cff(^) decreases monotonically 
with k, as shown in Fig. ||, this structure factor shows the 
fastest growth at the smallest wave number fc m ; n = 27r/L, 
allowed in a box of length L, and does not explain the 
dynamics of cluster growth. 



Next we consider the full theory of Sec. Ill with 
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Langevin noise included. Then the structure factor 
S„ n (k,t) — n 2 S^ n (k,r) is given by ( ^2|) as a sum over 
all pairs (A/u) of hydrodynamic modes, which exclude 
shear modes. Each term in this sum contains a factor 
exp[(£,\ + Cji) T ]> which leads to exponential growth when 
Ca + Cfi > 0. Inspection of the growth rates in Fig. || 
shows that this happens for (Xn) = (HH) for all k < k* H 
[defined in App. |a] below flA6)], and for (\{i) = (H+) 
for k below a certain threshold. The label A = + refers 
to the mode with vanishing rate constant as k — > (see 
App.g). 

To obtain a qualitative understanding of the cluster- 
ing instability, we present an analysis based on the most 
unstable pair of modes, (Au) = (HH). In this approxi- 
mation we obtain from (|22|), 



5+ n (fc,r) = n 2 WH n (y)wHn{-^)BHH{k) x 
exp[2Cg(fc)T] - 1 
2Cff(fc) 



(29) 



We analyze this expression in the dissipative range, 
klo <C 70, where Shh(^) is given in (p7), and where 
C,ji(k) ~ 70(1 — £ 2 fc 2 ). The component «)ff n (k) in ( [A^ ) 
of App. [a| vanishes for small k. Explicit calculation to 
the next order in k, as given in Ref. Jl4| ), shows that 
WHn(k) = —iklo/jo. These results combined with ( |l6| ) 
yield for the structure factor in the dissipative k range, 



S nn (k,t) = n 2 TxT + 



ikHl /exp[2 7o (l-e 2 fc 2 )r] -1 



27^ 



1 - ttk 2 



(30) 



This simple analytic approximation agrees within the 
small- A: range for which it is derived (say up to ka ^ 0.1 
m Fig. |), with the numerical solutions of the theory 
either with or without Langevin noise. Moreover, it 
demonstrates that the instability is driven through a cou- 
pling to the unstable 'heat' mode, which is, in the small- 
k range, a longitudinal velocity mode. The coupling of 
the density structure factor to the unstable heat mode is 
rather weak, 0(k 2 ), which explains why structure in the 
flow field appears long before density clusters appear. 

The wave number fc max (t) of the maximum growth of 
S nn in (|3^) determines the typical length scale of the den- 
sity clusters. For 270T 3> 1, it can be determined analyt- 
ically as L c \(t) ~ 27r/fc max (£) = 27r£|| v / 2tot, which is the 
same length scale, £||(t), as appeared in S\\. The good 
agreement between theory and MD simulations, shown 
in Fig. |^, confirm that the initial growth of density inho- 
mogencitics is indeed controlled by the longitudinal flow 
field with a length scale L c \(t) ~ J r/e at small inelas- 
ticity e = 2c?7o, and not by the transverse flow field with 
a length scale L v (t) ~ 2ir£± \J2^qt ~ y/r, independent of 
e. The pattern selection mechanism for the vortex struc- 
tures is very different from the mechanism that leads to 
the formation of density clusters. This is the more com- 
mon linear instability in density or composition fluctua- 
tions, which also occurs in spinodal decomposition H. 



V. SPATIAL CORRELATION FUNCTIONS 



Once the structure factors have been obtained, the cor- 
relation functions can be calculated by Fourier inversion. 
When a, b refer to n and T the components of S'oh(k, t) 
and G a b(r, t) are scalar isotropic functions only depend- 
ing on |k| and |r| respectively. When (a, b) = (a, (3) re- 
fer to Cartesian components u a of the flow field, then 
5 a /3(k, t) is a second rank isotropic tensor field, which 
can be separated into two independent isotropic scalar 
functions: 

S a0 (k, t) = k a kpS\\ (k, t) + (6 a/ 3 - k a kp)S±(k, t), (31) 

where S\\(k,t) and S±(k,t) are given by (||) with Sa and 
5b equal to ui and u± respectively. A similar separation 
applies to the velocity correlation functions, 

G a /3 (r, i) - V- 1 exp(zk • r)S Q/3 (k, t) 

k 

= r a fpG\\(r,t) + (S a/3 - r a rp)G±(r,t). (32) 



Firstly we note that the Fourier series in (32) can be re- 
placed by a Fourier integral, provided that the system 
is sufficiently large. Then, for periodic boundary con- 
ditions, as used in MD simulations, V^ 1 ^^ can be re- 
placed by (27r) -d / dk. Strictly speaking, isotropic sym- 
metry and separation of second rank tensor functions 
into two scalar functions only hold in the thermodynamic 
limit. 

Secondly, the inverse Fourier transform G(r, t) of 
S(k, t) only exists as a classical function, if S°°(t) = 
limfc^oo S(k, t) vanishes. If S(k, t) approaches a nonvan- 
ishing constant S°°(t), it yields a distribution <5(r). So, 



G(r,t) = S°°,5(r) 



dk 



(27T) 



z exp(ik-r)S + (k,t). (33) 



Note that the limiting values S%%, when expressed in 
terms of the rescaled variables, are given by the time 
and (wave number) independent values S a b(k,0). This is 
the reason for using the same notation as in Eq. ( |l6| ) . In 
the sequel it is convenient to also use the notation, 



G+(r,t) =G(r,i)-S°°(t)<S(r) 



(34) 



The structure factors S + (k, r) in (|2l|), calculated in the 
theory with Langevin noise, as well as S(k, r) in ( j24|) 
without noise are vanishing for large fc, and can be Fourier 
inverted. However, the functions S(k, r) in ( |l6| ) contain 
a part S(k, 0), which is independent of k in the relevant 
k interval, and which yields after Fourier inversion a con- 
tribution proportional to 6(r). These 'large k' contribu- 
tions are in fact the correlation functions of an elastic 
hard sphere (EHS) fluid as k — * 0, i.e. 
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G nn {r,t) ~n 2 T XT 6(r 
2T 2 (t) 



Gt(r,t) = 



G TT (r,t) 



dn 



S(r) 



G af3 (r,t) ~ !®5(r)6 af ,. 



(35) 



According to Sec. El 'large k' here means y^fo/lo <C 
k <C min{27r/io, 27r/cr}. Here G nn is the coarse grained 
density-density correlation function for EHS, in which 
the Fourier components with ka > 2tt have been dis- 
carded. In App. [B]we derive the formulas, necessary for 
the analytic and numerical Fourier inversion of S + (k, t), 
as defined in (1341). 



A. Incompressible limit 



There is an interesting limiting case of the theory, the 
incompressible limit 0] , that greatly simplifies and eluci- 
dates the analytic solution of the full set of coupled lin- 
earized equations (0) for hydrodynamic fluctuations. It 
is well known from fluid dynamics and the theory of tur- 
bulence j32 19 that ordinary elastic fluid flows are quite 
incompressible, which implies that V • u = and that 
the longitudinal mode u/(k, t) ~ 0. Then, the nonlin- 
ear Eq. (|^) for the transverse flow field or, equivalently, 
for the vorticity, practically decouples from the remain- 
ing hydrodynamic equations. In the comoving reference 
frame there is only a nonlinear coupling of the tempera- 
ture fluctuations to the transverse flow field through the 
nonlinear viscous heating, ?y|Vu| 2 . We therefore expect 
that the IHS fluid in the nearly elastic case can be consid- 
ered as incompressible, at least to lowest approximation. 

What are the consequences of this assumption for the 
Fourier modes? The structure factor S nn (k,t) of den- 
sity fluctuations does not evolve in time. The tempera- 
ture fluctuation 5T(k, t) in (Q) simply decays as a kinetic 
mode and the average temperature stays spatially homo- 
geneous. Clearly, the assumption is too drastic a simplifi- 
cation to describe the clustering instability. However, an 
approximate theory based on vorticity fluctuations alone 
is justified to de scri be the patterns in the flow field, as 
discussed in Sec. [iy| 

So, we combine the assumption of incompressibility, 
S\\{k,t) = 0, with S±(k,t) in @, using @ and (§§). 
This enables us, for thermodynamically large systems, to 
explicitly calculate the correlation functions G a p(r, t) of 
the velocity field, by inverse Fourier transformation, i.e. 



Ga/sM) = f a rpG+{r,t) + (5 af} - f a f p )G+(r 7 t) 



dk 



— — t exp(ik • r)(S af3 - k a k )SX(k, t). (36) 

(Z7TJ 



The Fourier transform has been calculated in App. 
and yields for the two scalar functions G\ (r, t) with 
A = {||,T} 



T(t) 



(37) 



where g\(x,s) is given in (Bll) of App. EL and 70 and 
£1 are defined below Eqs. (|) and (0) respectively. 

As an explicit example wc show the result for inelastic 
hard disks, which is most relevant for a comparison with 
existing computer simulations, in Fig. ||. Their analytic 
form is, 

1 f s 

5||(z,s) = 7^ J ds'exp(s')[l-exp(-a; 2 /4s')] 

g±(x,s) = -g\\(x,s) + J ds cxp(s ) — . (38) 

The function g±(x,s) has a negative minimum, while 
gu (x, s) is positive for all x, s, d; there are algebraic tails 
gii [x, s) ~ — (d — l)g±(x, s) ~ x~ d with a correction term 
of C(exp(— x 2 /4s)), explicitly given in (B12). Similar 
algebraic tails occur in noncquilibrium stationary states 
in driven diffusive systems ]33|,|34j]. These functions have 
structure on hydrodynamic space and time scales where 
both x = r/£j_ and s = 2jqt can be either large or small 
with respect to unity. At small inelasticity (70 — > 0) the 
dynamic correlation length £j_ and mean free path Iq are 
well separated. 

A more systematic comparison between the theoreti- 
cal predictions (Bq) and molecular dynamics simulations 
is made in Ref. In Fig. |^, we show the results from 

a single simulation run at packing fraction <f> = 0.4, and 
small inelasticity a — 0.9. The parallel part G\\(r,t) ex- 
hibits a tail - r~ d (see Fig. 2(a) in Ref. 0) and shows 
good agreement, well beyond the crossover time r c (de- 
fined in Fig. [2]) that separates the linear regime from the 
nonlinear clustering regime. The minimum in G±(r, t) at 
L v (t) can be identified as the mean diameter of vortices, 
shown in Fig. |l]. The analytic result for G±(r, t) in (|3"§| ) 
for large times shows that that L v (t) ~ 27t£_i_- v /27ot is 
growing through vorticity diffusion. 

Apart from the restrictions to hydrodynamic space and 
time scales, there are two essential criteria limiting the 
validity of the incompressible theory: (i) System sizes 
L must be thermodynamically large (L ^> 2tt^±), so that 
Fourier sums over k-space can be replaced by k-integrals. 
(ii) Times must be restricted to the linear hydrodynamic 
regime (t < r c ), so that the system remains close to the 
HCS. It appears that our description of the fluctuations 
in terms of a Langevin equation based on incompressibil- 
ity is confirmed by the simulations in the linear regime 
r < r c for small inelasticities. 



B. Compressible Flows 

In this subsection we extend the theory to compress- 
ible flows H. The description of the velocity fluctu- 
ations G a p{Y,t) in the previous subsection was based 
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on fluctuating hydrodynamics for the vorticity fluctua- 
tions only, i.e. the absence of longitudinal fluctuations 
(incompressibility assumption). Fig. ||(b) confirms that 
this assumption is very reasonable indeed, as S^(k,t) — 
S\\(k,t) — T/mn is vanishingly small down to very small 
k values. However for the smallest wave numbers, the 
incompressibility assumption breaks down. As the anal- 
ysis of ( p6|) and (E8|), as well as the numerical eval- 
uation in Fig. [3] show, the structure factor S\\ (k —> 
0, t) = S±(k — *■ 0,i). This implies for large distances 
G a p(r,t) ~ S±(k — > 0,t)6 a p5(r), and thus the absence 
of algebraic long range correlations on the largest scales 
(r 3> 27r£||). Therefore, we can already conclude that 
the asymptotic behavior of G±(r,t) and Gn(r,t) cannot 



be r 



Instead the 



tail, obtained in the previous 



subsection, describes intermediate behavior which is ex- 
ponentially cut off at a distance determined by the width 
of St(k,t). This width can be estimated from the eigen- 
values of the hydrodynamic matrix, more precisely from 
the dispersion relation of the heat mode, which is a pure 
longitudinal velocity ui for k — > 0. To second order in k 
its dispersion relation is given by Cn{k) = 70(1 — fc 2 £jj), 
where £|| is given in Eq. (|A7|) of App. |aJ Note that 
for small inelasticities £|| and £j_ are well separated, as 

£|| ~ V e > whereas £1 ~ 6 ~ 6r ~ l/v^- 

Using approximation (|2^) for S\\(k, t), the structure 
factor S^p(k,t) can be written as 

S+(k,t)«^ / d S 'exp( S ') \k a kpexp(- S 'k 2 ^) 



+(6 a p - k a kp) exp(-s'/c 2 ^) , 



(39) 




where s — 2j t. If the system is thermodynamically 
large (L ^> 27r£||), Gt(r,t) and G\_{r,t) can be obtained 
by performing integrals over k space and are expressed 
as integrals over simple functions, as derived in App. |^. 
Here we only quote the results for d = 2: 



(40) 



for A =||, -L, where x\ — r/£\, eni = 1 and a± = —1. We 
first observe that, in the time regime t ~ I/70, G^(r,t) 
has structure both on the scale r ~ 2tt^± as well as on 
r ~ 27r^||. Moreover, Gjj" (r, t) is positive both in the in- 
compressible as well as in the compressible case because 
< In Fig. H we show the above approximations 
for different values of the ratio together with the 

incompressible limit result of the previous section, which 
is obtained for £|| — > 00. At finite e, Eq. po] ) describes 
exponentially decaying functions at distances r > 27r£ii. 
Moreover, upon increasing the inelasticity the minimum 
in G±(r, t) becomes less deep and vanishes if £|| =£j_. 



The predicted spatial velocity correlations G11 (r, t) and 
Gx (r, t) have been obtained by performing inverse Bessel 
transformations on the numerical results for S\\(k,t) and 



S±(k,t). At small inelasticity 



> 



0.9) the functions 



Gu(r, t) and G±(r,t), calculated from the full set of hy- 
drodynamic equations, differ for r < 27r£|| only slightly 
from the results for incompressible flow fields (see dis- 
cussion in the previous section). However, the algebraic 
tails ~ r~ d in Gn(r,t) and G±(r,t) for r > 2ir£±, as de- 
rived in the previous subsection, are exponentially cut off 
for r > 27r£||, as implied by Eq. ((40|). As the correlation 

lengths £j_ ~ l/v^ an d £|| ~ V e are wen separated for 
small e, there is an intermediate range of r values where 
the algebraic tail ~ r~ d in Gu (r, t) can be observed. 

At higher inelasticity £y and £j_ are not well separated 
and, as a consequence, there does not exist a spatial 
regime in which the longitudinal fluctuations in the flow 
field can be neglected and the regime of validity of the 
incompressible theory has shrunk to zero. Fig. [l^ com- 
pares results from incompressible and compressible fluc- 
tuating hydrodynamics with simulation data for G± (r, i) 
at <f» = 0.4 and a — 0.6, and confirms the necessity of 
including longitudinal velocity fluctuations to calculate 
the spatial velocity correlations at reasonably large in- 
elasticities. 

Fig. [ll] shows the spatial correlation functions G nn 
and G n T , which are the inverse Fourier transforms of the 
structure factors S a b(k,t), shown in Fig. fj]. The spatial 
density correlation G nn (r,t), obtained numerically from 
Snn(k, t), exhibits a negative correlation centered around 
a distance which for large times grows as y/r. Compari- 
son with simulation results confirms that the present the- 
ory correctly predicts the buildup of density correlations 
in the linear time regime r < t c . For a more comprehen- 
sive comparison of the velocity and density correlation 
functions with MD simulations at different densities and 
different inelasticities, we refer to Ref. 0, where also the 
range of validity of the present theories will be investi- 
gated. 

In summary, the typical length scales, of the vortices 
L v (t) ~ Jt and of the density clusters L c \{t) ~ \frji 
in Sec. |v|, also correspond to the typical length scales 
on which the equal-time correlation functions G± (r, t) , 
G\\{r,t) and G nn (r,t) show structure. In zeroth approx- 
imation the HCS is an adiabatically cooling equilibrium 
state with short range correlation functions, given in (|3^) . 
As time increases, long ran ge c orrelations develop which 
are at most of 0{e d ^ 2 ) [see ( p7| ) and (fiC|)1, and extend far 
beyond the microscopic and kinetic scales. 



VI. CONCLUSION 

In this paper the structure factors S , a b(k, t) and corre- 
sponding spatial correlation functions G a i,(r, t) have been 
calculated and compared with 2D molecular dynamics 
simulations for weakly inelastic hard disk systems, where 
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e = 1 — a 2 is small. For strongly inelastic systems the 
macroscopic equations are not known. For weakly inelas- 
tic systems we have assumed the hydrodynamic equa- 
tions of an elastic fluid, supplied with an energy sink 
representing the collisional dissipation. Also, we have 
assumed that the homogeneous cooling state (HCS) is 
an adiabatically cooling equilibrium state, which is o nly 
correct to lowest order in e. In fact, in Ref. [ [28|]30| 
the Enskog-Boltzmann equation is used to calculate e- 
dependent corrections to the velocity distribution in the 
HCS. As follows from the Enskog-Boltzmann equation, 
the pressure for inelastic hard spheres decreases linearly 
with decreasing a. However, in our lowest order the- 
ory it has been set equal to the pressure of elastic hard 
spheres (EHS). Moreover, the collision frequency u> and 
the transport coefficients of viscosity and heat conductiv- 
ity have been assumed to be given by the Enskog theory 
for EHS, but there are substantial corrections for inelas- 
tic hard spheres. Moreover, there appear new transport 
coefficients, which are absent in EHS fluids pGLp7|. 



The most important parts of this paper are Sees. Ill 
and |y|, where the basic theory is developed, modeled 
on the Cahn-Hilliard theory for spinodal decomposition. 
Based on the idea that the HCS for weakly inelastic hard 
spheres is essentially an adiabatic equilibrium state with 
a slowly changing temperature, we have formulated a 
fluctuation-dissipation theorem for this state, which en- 
abled us to construct Langevin equations for hydrody- 
namic fluctuations. The full set of coupled equations 
( |ll| ) for the structure factors can only be solved numer- 
ically. To understand the physical excitations that drive 
the instabilities, we have presented a theoretical analysis 
of the structure factors using a spectral analysis of the 
unstable Fourier modes. The structure factor S±(k, t) for 
the transverse flow field can be calculated exactly. For 
the longitudinal one, S\\ (k,t), we have obtained the sim- 
ple analytic approximation (|28| ) . It only accounts for the 
dominant contribution of the heat mode, and gives an 
almost quantitatively correct description of S\\(k, t) for 
all times, but it is missing the little dip which appears 
both in the numerical results of Fig. 0(a), as well as in 
the MD simulation results of Fig. ^(b). 

The dynamics of the transverse and longitudinal flow 
fields on the largest length scales are controlled by stable 
purely diffusive modes with very different diffusivities. 
There are no sound modes on these length scales. The 
structure factors S±(k,t) and Su(k,t) are bounded at all 
k and t by their initial values, and they develop spa- 
tial structure on length scales L v (t) ~ 2tt£±^/2jot and 
Z/||(i) ~ 27r£||- s /27 T respectively. 

Agreement between the predictions of fluctuating hy- 
drodynamics with Langevin noise for iS^ and S±, and the 
results of MD simulations is very good. 

Calculation of the structure factor S nn (k,t) for the 
density fluctuations is essentially a linear stability ana- 
lysis, which describes the early stages of clustering. It 
is expected to break down, because the density fluctua- 
tions are predicted to grow at an exponential rate. This 



is indeed seen to happen as the time r approaches the 
crossover time r c to the nonlinear clustering regime, as 
defined in the caption of Fig. ||. The agreement with MD 
simulations for r < r c is again quite good. Moreover, 
the simple analytic long wavelength approximation ( |30| ) 
agrees in the most relevant small-fc range almost quan- 
titatively with the numerical results of the theory cither 
with or without Langevin noise. 

One would expect that the agreement between the- 
ory and simulations would also break down for Su and 
Sj_ when r approaches r c . However, it has been shown 
that 5|| and S±, calculated from this linear the- 
ory, remain in agreement with the MD simulations for 
r > t c . More surprisingly, the energy decay Ejt) in Fig. 
H has been quantitatively explained in Ref. fllif over the 
whole simulation interval (r < 150 ~ 2r c ), using struc- 
ture factors S± and Su , calculated from the present linear 
theory. 

The density structure factor S nn (k,t) shows that spa- 
tial density fluctuations in undriven IHS fluids are unsta- 
ble, and lead to the formation of density clusters. The 
linear instability is driven by longitudinal velocity fluc- 
tuations and described by a coupling coefficient of 0(k 2 ) 
in S nn . The fluctuations in the flow field are only rel- 
atively unstable, and do not lead to exponential growth 
of the corresponding structure factors. Nevertheless, the 
dissipativc IHS fluid develops structure on intermediate 
scales with typical length scales L c \(t) ~ £,\\V^t ~ \J T I ( - 
for the mean cluster sizes, and L v (t) ~ S,±^/er ~ \fr for 
the mean vortex diameters. 

In the literature the clustering instability has also 
been studied on the basis of the nonlinear hydrodynamic 
Eqs. which show that clustering in the late stages of 
evolution is driven by the viscous heating term, ?y|Vu| 2 . 
The typical length scale of the clusters is then related to 
the shear mode, £j_ = y v /ojjq , where v is the kinematic 
viscosity and oj the collision frequency. 

Section [v| deals with spatial correlations. The assump- 
tion of incompressible flow for nearly elastic hard spheres 
(a > 0.9) — i.e. based on the relative instability of shear 
modes only — leads to spatial velocity correlations, in- 
cluding algebraic r~ d tails, that are correct up to large 
distances (r < 27r£ii). We have shown by explicit cal- 
culation that at small inelasticities S^(k,t) essentially 
vanishes for all wave numbers except at very small k val- 
ues (k < 1/£||)> where the assumption of incompressible 
u fluctuations, made in subsection V.A, breaks down. 
Consequently, at small inelasticities the most important 
qualitative modification that adds to the spatial cor- 
relation function G\\ (r, t) of subsection V.A is to provide 
an exponential cutoff for the r~ d tail at the largest scales 
r > 2-7r£j|. At larger inelasticities the nonvanishing con- 
tributions from St (k, t) modify G\\ (r, t) and Gj_(r, t) sig- 
nificantly at all distances. 

The good quantitative correspondence between the- 
ory and computer simulations shows that our theory for 
structure factors, S a p(k, t) and S nn (k, t), and spatial cor- 
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relation functions, G a p(T,t) and G nn (r,t), is correct for 
wave number, position and time dependence in the rele- 
vant hydrodynamic range and for inelasticities (a > 0.6) 
that are not too large. 

Moreover, we have emphasized that there exist two dif- 
ferent mechanisms for structure formation and pattern 
selection, active in granular fluids. The first one, driv- 
ing the formation of vortex structures, is the mechanism 
of noise reduction [14j], comparable to peneplanation in 
structural geology with Ayers Rock (Mount Uluru) as a 
spectacular example. This mechanism selectively sup- 
presses velocity fluctuations at shorter wavelengths as 
compared to longer ones, as well as longitudinal fluc- 
tuations as compared to transverse ones. The second 
mechanism, driving density clustering, is the more com- 
mon instability in density or composition fluctuations as 
occurring in spinodal decomposition. In freely evolving 
granular fluids the density is most unstable at typical 
wavelengths L c \(t) ~ 2-7r/fc max (£). 
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APPENDIX A: APPENDIX: TRANSPORT 
COEFFICIENTS 

Linearization of the hydrodynamic equations (|^) 
around the HCS results in the following set of equations: 



dtSn — —n'V ■ u 

d t u = --Vp + vV 2 u + ivi - v) 
P 

dtS T=^V 2 ST-^-V-u-6T. 

an an 



VV u 



(Al) 



The pressure p is assumed to be that of elastic hard 
spheres (EHS), p = nT(l + Vt d xna d /2d), where fl d = 
27r rf / 2 /r(d/2) is the d-dimensional solid angle, and x( n ) 
the equilibrium value of the pair correlation function of 
EHS of diameter a and mass m at contact. The kine- 
matic and longitudinal viscosities v and vi, as well as 
the heat conductivity k are also assumed to be approx- 
imately equal to the corresponding quantities for EHS, 
as calculated from the Enskog theory |24|] , where pv = r\ 
and pv\ — 2n(d—l)/d + ( are expressed in shear viscosity 



r\ and bulk viscosity £. The collisional energy loss in ([ 
r = 2^qujT, is proportional to the collision frequency 



irm 



(A2) 



as obtained from the Enskog theory. In three dimen- 
sions we use the Carnahan-Starling approximation \ — 
(2 — <f)/2(l — cf)) 3 p5| , in two dimensions the Verlet- 
Levesque approximation ^ = (1 — 70/16)/(l — 0) 2 [j36f . 

In the body of the paper we need the eigenvalues Ca (k) 
of the asymmetric matrix M a j, (k) , introduced in (Js[) , and 
its right and left eigenvectors, which are obtained from 



M(k)w A (k) =CA(k)w A (k) 
M T (k)v A (k) =CA(k)v A (k), 



(A3) 



where M T is the transpose of M. Here A = ± labels the 
sound modes, A = H the heat mode, and A =_L labels 
(d — 1) degenerate shear or transverse velocity modes. 
The eigenvectors form a complete biorthonormal basis, 
i.e. 

^w Aa (k)w Aia (k) = (v A |w M ) = <5 Am 

a 

£|w A (k))(v A (k)|=I, (A4) 



where (.|.) is the usual scalar product. The relations ( A4) 
allow a spectral decomposition of M(k) in the form 



M ab (k) = ^% a (k)CA(A;)UA6(k). 



(A5) 



Moreover the eigenvalue equation, det[£(fc)I— M(k)] = 0, 
is an even function of k. Consequently, M(k) and M(— k) 
have the same eigenvalues, which are either real or form a 
complex conjugate pair. So we choose CaW = Ca(~ k) = 
C\(k). The corresponding eigenvectors of M(— k) in 
case of propagating sound modes are obtained from the 
transformation {w + (— k), v_ (— k)} <-> {w_(k),v + (k)}. 
All other eigenvectors, corresponding to nonpropagating 
modes, are invariant under the transformation k — * — k. 

The right and left eigenvectors of the hydrodynamic 
matrix at zero wave number are given by 



Wj.(k) 

w ff (k) 
w + (k) 
w_(k) 



(0,0,0,1) 
(0,0,1,0) 
(l,-g(n),0,0) 
(0,1,0,0) 



v x (k) = (0,0,0,1) 
v H (k) = (0,0,1,0) 
v+(k) = (1,0,0,0) 
v_(k) = (g(n), 1,0,0). 



(A6) 



They characterize the (d — 1) shear modes (A =-L), the 
heat mode (A = H) and the two sound modes (A = ±). 
Higher- k corrections are given in Ref. |14|j . The shear 
or transverse velocity mode is decoupled from the other 
modes, and has a growth rate Cl(^) = 7o(l — ^ 2 Cl)- 

In the dissipative range (fc^o -C 7o), as k — > 0, the heat 
mode is a pure longitudinal velocity fluctuation, while the 
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sound modes are a mixture of density and temperature 
fluctuations. To first order in k, density and tempera- 
ture fluctuations couple to the heat mode, and longitudi- 
nal velocity fluctuations couple to the sound modes. For 
large wave numbers (klo S> 70) the conventional char- 
acter of heat and sound modes is recovered. Besides 
the (uncoupled) shear mode, the heat mode is unsta- 
ble for k < k* H , where k^ is the root of Cff(^) = 0; or 
equivalently the root of det[M(k)] = 0. This yields the 
stability threshold, k* H = [g{n)xTP - 1]^ 2 /^t ~ h/y/e- 
Its dispersion relation to second order in k is given by 
( H (k) =70(1-^), where 



I 2 



1 



tiTxt 



( V \ L , nd\ P_ 

V nT J \ x 9n dnT 



(A7) 



Note that £11 diverges as 1/e, while £j_ ~ l/v^i as a 
consequence the correlation lengths £|| and £j_ are well 
separated for small inelasticity ||]. 

The sound modes are stable for all wave numbers. In 
the dissipative regime (klo "C 70), they correspond to 
nonpropagating modes. The mode labeled A = + is for 
small k a linear combination of a density and a temper- 
ature fluctuation with (+(k) ~ — 7ofc 2 £ 2 _, and the mode 
labeled A — — is a pure temperature fluctuation, with 
(-(k) ~ 7o(— 1 + & 2 £ 2 ), which is a kinetic mode. To 
lowest order in k the new correlation lengths are given 
by 



2 (nT7 



1 4. H^K 
X 9n 



1 



e(k) = 



2 7 2 KnTJ 



tiTxt 

1 + H^K ' p 
X dn 



dnT 



(A8) 



At a wave number of the order of Jo/Iq, their dispersion 
relations merge and the modes start to propagate. 



APPENDIX B: APPENDIX: FOURIER 
TRANSFORMS 

To calculate the tensor velocity correlation function 
G Q( a(r, t) by Fourier inversion from S a p(k,t), we start 
from Eqs. @, ©, © and @, and consider first the 
incompressible limit where S'||(fc,i) = 0, i.e. 



dk 



(27T) 



^ exp(zk • r)(6 a p - k a kp)Sj_(k, t). 



(Bl) 



According to (|3^) G a p can be split into two scalar func- 
tions, G|i and Gj_, which will be expressed in Sj_. The 
simplest functions to calculate are the trace and parallel 
part of G a 0, i.e. 



Gp P (r, t) — G+ Q (r, t) 



dk 



G+(r,t)=f a fpG+ (v,t) 



dk 



^-pexp(zk.r)[l-(k.f)l5]:(fc,i). 



(B2) 

To carry out the ci-dimensional angular integrations for 
d > 2 we express the infinitesimal solid angle as 



dk= (sinfli) 



d-2 



(sin0 rf _ 2 )d0i...d0 rf _ 2 d0, (B3) 



where 8 n € (0,7r) are polar angles and (f> £ (0,27r) is an 
azimuthal angle, and we note that the full solid angle is 
Q d = / dk = 27r rf / 2 /r(d/2). Then we use the relation 

dk C d0(sin 9) d - 2 exp(i/cr cos 6) 

— exp(>k • r) = -y . , 

S2d J d6>(sm6>) d 2 

2 \ d/2-1 

r(d/2)J d/2 _ 1 (fcr), (B4) 

where the integral representation (8.411.7) of Ref. p7| ] 
has been used for the Bessel function J v (z). Then Eqs. 
(B2) become 



Gp P (r, t) — 
G+(r,t) = 



1 



" (2 7r )<i/2 r d/2-i y Q 
d- 1 



( 27r )d/2 r d/2 J Q 



dfcfe d / 2 J d/2 _ 1 (fcr)S+(A;,t) 

dfcfc d / 2 - 1 J d/2 (fcr)^+(fc,t). 

(B5) 



With the help of the recursion formula for Bessel func- 
tions, zdJ v (z) / &z + vJ v (z) = zJ v -x(z), together with the 
general relation 

G+(r,f) = [G+(r,t) - G+(r,t)]/(d- 1), (B6) 
we obtain G*(r, i) from Gjl" (r, t), i.e. 

G+(r,<) = G+(r,i) + f^-) ^M). (B7) 

This is a well known relation in the theory of homoge- 
neous and isotropic turbulence in incompressible flows 
(see Refs. || and g§, Chap. 3). 

In the general case Su (k, t) is nonvanishing and we have 
from ( |3l] ) an additional part, denoted by G a p(r, t), com- 
ing from k a kf}S«(k, t). Here we have the relations 



G pp (r, t) 
G + Ar,t) 



dk 

W) 

1 f dk 



1 



1 exp(ik-r)Sj{k,t) 



1 exp( l k-r)[l~(k-r) 2 }S+(k,t). 

(B8) 



(27T) 
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The results for these functions can be read off directly 
from (B2) and (B5), and the parallel part is obtained 
from (IBST) as 



Gj(r,t)=G; p (r,t)-(d-l)Gl(r,t) 



In the incompressible limit, where S\\(k,t) — 0, the ve- 
locity correlations in the fl ow field can be cal cula ted by 
Fouri er inversion of ( B13| ). Comparison of (B12) with 
(|B10|) shows that 



Gl(r,t)+r—G+(r,t). 



(B9) 



(B14) 



The derivation of the second line parallels that of Eq. 
(B7). Note also that the role of Gu and G± has been 
interchanged with respect to the incompressible case. 

Fourier inversion of any of the scalar functions S a b(k, t) 
with a,b = \n,T} is covered by the first line of Eqs. 



where 



/ \ 9g\\(x,s) d-1 (d x 2 



( |B2[) and (|B5| ). We consider first G\\ and G± in the in- gx(x,s) = exp(s) 



compressible limit, where S±(k,t) is given by (|26|). The 
Cahn-Hilliard theory in the incompressible limit will be 
considered later. 

Inspection of ( p6|) shows that the large-fc limit of 
Sj_(M) is S^ 1 = T(t)/mn, leaving S±(k,t) as a remain- 
der. This may be written as 



cxp(— £ 2 /4s) 

(47Ts) d / 2 



1 



j 9ll (x,s). (B15) 



S+(k,t) 



T(t) 



ds'exp[(l-fc 2 el)s']. (BIO) 



Using Eq. (B5) for the parallel velocity correlation func- 
tion Gf(r,t) and (B7) to determine C\_(r,t), we obtain 

G+(r,t) = (T(t)/mnti)g x (x,s) for A = {||,_L} with 
s = 2"fo t and x = r/£j_, which is is valid in dimensions 
d > 2. Moreover, g\(x, s) is given by 



d-l 



9\\( x > s ) = ZdlZd L ds'exp( S ')7^ 2 , 4s , 



d x 2 



9±(x,s) 



2 n d/2 x d J q 

A / / ^exp(-x 2 /As') 

^ -Ms) (w)d/2 



d-V 



9\\(x,s). 



(Bll) 



The Bessel transform (B5) of exp(— /3k ) in (Bid) has 
been calculated using Eq. (6.631.5) of Ref. |37f1 , where 
j(a,z) — J Q dt exp(— t)t Q_1 is the incomplete gamma 
function. For d = 2 it reduces to 7(1, z) = 1 — exp(— z) 
and for d — 3 to 7(3/2, z) — ^/n(j)(y/z)/2 — y / iexp(— z), 
where cj>(z) is the error function. We observe that gn (x, s) 
is positive for all x, s. For large distance, x 2 3> 4s, the 
functions g\(x, s) show algebraic tails ~ l/x d . This can 
be seen by noting that j(a,x 2 /4s) approaches T(a), so 
that 

S||(a;_L.,s) =-{d- l)g±(x ± ,s) - (^Jr^J I ex P( s ) _ x ]> 

(B12) 



where Qd is defined below (B3). 

The theory without noise for the transverse structure 
function is according to the discussion in Sec. IV given 
by 

Sx(M) = St(k,t) = I^exp[2 7 o(l - fc 2 Ci)r]. (B13) 
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FIGURE CAPTIONS 



1. Left: Velocity field after r = 80 collisions per particle. The density is then still more or less homogeneous. Right: 
Density field at r = 160. System of inelastic hard disks at a packing fraction <fi = \na 2 = 0.4 and a = 0.9. 

2. Kinetic energy per particle E versus number of collisions per particle r for ip = 0.4 and a = 0.9. Initially E is 
equal to the temperature T and follows Haff 's homogeneous cooling law (^j). The arrow indicates a crossover time 
t c ~ 70 from the homogeneous cooling state to the nonlinear clustering regime. Then spatial inhomogeneities 
become important and slow down the energy decay. The dashed line represents Haff 's law (^) and the dot-dashed 
line the result of Ref. [|l4| for the long time energy decay. 

3. Growth rates Ca/7o for shear (A =-L), heat (A = H) and sound (A = ±) modes versus ka for inelastic hard disks 
with a = 0.9 in (a) and a = 0.6 in (b) at a packing fraction cf> — 0.4 (Iq ~ 0.34ct). The shear and heat mode 
are unstable for k < k\_ and k < k* H respectively. Imaginary parts of the sound modes, indicated by the dashed 
lines, vanish for Hq <C 70 ■ 

4. Ratio versus packing fraction 4> of inelastic disks, with definition of £|| in (A7) and of £j_ below (0). 

5. Density structure factor S nn , in units l/a 2 , versus ka for <j) = 0.4 and a = 0.9, at t = 10,20,30 and 40 
collisions per particle, exhibits the clustering instability with a growing maximum at k max (t), which shifts to 
the left. Solid and dashed lines are the numerical solutions of (Q) and (|24|) with and without Langevin noise 
respectively. They differ appreciably, except at small k. The simple analytic approximation (|3p|), shown in (a) 
as dot-dashed lines for r = 30 and 40, gives a good description in the long time and long wavelength limit, (b) 
Numerical solution of ( pr| ) compared with MD simulation results (courtesy J.A.G. Orza et al. (J]). 

6. Structure factors of velocity fluctuations S\\ and S±, in units Tqct 2 /to, versus ka for 4> — 0.4 and a = 0.9 
illustrate the phenomenon of noise reduction at small wavelengths, (a) 5y at r = 10, 20, 30 and 40, where 
solid lines represent the numerical solution of (|l4|) , and the dot-dashed lines represent the approximate analytic 
result (ps|). The dashed line represents the numerical solution of the 'noiseless' Cahn-Hilliard theory (|24|) at 
t = 10, which deviates substantially from the solid line at r = 10, except near k = 0. (b) Comparison with MD 
simulation results (courtesy of J.A.G. Orza et al. pj) at r = 20 for S± (squares) and S\\ (circles). The numerical 
solutions of (H) and ((2J) with (solid lines) and without Langevin noise (dashed lines) approach each other for 
long times and long wavelengths. 

7. Rescaled structure factors (a) Stt, (b) S n r, (c) Im(S n i), and (d) Iui(Sti), all in units a 2 , versus ka for same 
parameters as used in Fig. g. Both S n i and Sti change sign under k — ► — k and are therefore purely imaginary. 
The structure factors in (b), (c) and (d) vanish initially and develop structure as time increases. Stt develops 
structure on top of its initial (plateau) value 2/dn. 

8. Comparison of theoretical predictions (|38|), based on incompressibility, with MD simulation results for Gn (filled 
circles) and G± (open circles), in units To/m, versus r/a At <j) = 0.4 and a — 0.9. 

9. (a) g\\(x, s) and (b) g±(x, s) versus x = x± — r/£j_ for s = 2joT = 2. The solid lines corresponds to (|38|) in 
the incompressible limit - * °°), and the dashed lines to approximations (f4(i|), for £||/£± — lj 2, 5, 10. As 
£||/£-L decreases, the r~ 2 tail in (a) is cut off exponentially at smaller distances and finally disappears at £|| = 
The depth of the minimum in (b) decreases with decreasing anc ^ finally disappears at £|| = 

10. Perpendicular velocity correlation G±, in units To/m, versus r/a for packing fraction <fi = 0.4, relatively high 
inelasticity a — 0.6 and r = 40. Simulation results are compared with prediction ( pq ) of the incompressible 
theory (dashed line), and the numerical solution (solid line) of the full set of fluctuating hydrodynamic equations. 

11. Spatial correlation functions (a) G nn (r,t), in units 10~ 3 /er 4 , and (b) G n T{r,t)/T(t), in units 10~ 3 /er 2 , versus 
r/a obtained numerically from the structure factors shown in Fig. at the same parameters as used in Fig. ||. 
Both functions show a growing correlation length. 
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